
% comparison of heterogeneous and homogeneous hyps learners:
% check how well each does over a range of hyperparameters

addpath('gpml-matlab-v3.1-2010-09-27');
startup;

n_expts = 10;
n_reps = 4;

gendata(); % so n is correct, mostly
func = @covSEiso;
D = load('gp_data'); % much of this is CVed from cross_validate.m
dim = size(D);
n = size(D,1);
% ntot=n;
%D = D(1:ntot,:); % throw away data if need be
dim = dim(2)-1; % i.e. data are from a dim-dimensional space; the -1 is not counting the class label
n_class = max(D(:,dim+1));
K = zeros(n,n,n_class);sigma_noise = 1e-7;

hyps = load('gp_hyps')'; hyps = reshape(hyps, 1, n_class*dim);
approxF = zeros(n*n_class);
l = -30;

results = zeros(n_expts, 4); % heterogeneous & error; homogeneous & errors
data = zeros(n_reps,2); % how well hetero does each time; how well homo does
for expt = 1:n_expts
  fprintf('\nExpt: %i\n', expt);
  for trial = 1:n_reps
    fprintf('Trial: %i\n',trial);
    gendata([0; 4 * (0.5 * (n_expts + 1) - expt) / n_expts]);
    D = load('gp_data'); % much of this is CVed from cross_validate.m
    X = D(:,1:dim);
    y1 = D(:,dim+1);
    y = reshape(kron(ones(n,1),1:n_class),n*n_class,1)==repmat(y1,n_class,1); % expressed as 1-of-n_class encoding
  
    hyps = minimize(hyps, @eq_3_44, l, func, n, n_class, X, y, approxF);
    dlmwrite('gp_hyps', hyps); % fingers crossed
    disp('Step 1')
    cross_validate()
    tmp = load('cv_results');
    data(trial, 1) = tmp(1);

    Hyps = reshape(hyps, str2num(func()), n_class);
    Hyps = sum(Hyps')'/n_class;
    disp('Step 2')
    Hyps = minimize(Hyps, @eq_3_44, l, func, n, n_class, X, y, approxF);
    dlmwrite('gp_hyps', repmat(Hyps,n_class,1)); % fingers crossed

    disp('Step 3')
    cross_validate()
    disp('Step 4')
    tmp = load('cv_results');
    data(trial, 2) = tmp(1);
  end
  av = sum(data(:,1)) / n_reps;
  results(expt,1) = av;
  vars = (data(:,1) - av) .^ 2;
  results(expt,2) = sqrt(sum(vars)/n_reps) / sqrt(n_reps);
  
  av = sum(data(:,2)) / n_reps;
  results(expt,3) = av;
  vars = (data(:,2) - av) .^ 2;
  results(expt,4) = sqrt(sum(vars)/n_reps) / sqrt(n_reps);
  dlmwrite('expt3_results', results);
end
results = results/n; % so it's all fractional
% I'm pretty sure it's simple division like that
dlmwrite('expt3_results', results);
disp('Done');